#pragma once
#include <cuda_runtime.h>
#include <glm/glm.hpp>

namespace PhysLeo {

/**
 * static peridynamic nonlinear material.
 * \f$ g=h=\sum^4_{i=1}{p_iA_i} , G=1\f$.
 * \f$A_i\f$ is A class nonlinear basis function, \f$p_i\f$ is the nonlinear parameters which controls the combination of different nonlinear basis.
 */
template<typename T>
class PdmNonlinear {
public:
    /**
    * g function in Pdm nonlinear material is \f$ \sum^4_{i=1}{p_iA_i} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T gFunction(T x, T* ptr_material)
    {
        T nonlinear_params[4] = { ptr_material[2],ptr_material[3] ,ptr_material[4] ,ptr_material[5] };

        auto x2 = x * x;
        auto x3 = x2 * x;
        auto x4 = x2 * x2;
        auto x5 = x4 * x;
        auto sum = nonlinear_params[0] * (0.5f *x2 - log(x) - 0.5f);
        sum += 0.5f*nonlinear_params[1] * (1.0f / 3.0f*x3 + 1.0f / x - 4.0f / 3.0f);
        sum += 1.0f / 3.0f*nonlinear_params[2]*(1.0f / 4.0f*x4 + 1.0f / 2.0f / x2 - 3.0f / 4.0f);
        sum += 1.0f / 4.0f*nonlinear_params[3]*(1.0f / 5.0f*x5 + 1.0f / 3.0f / x3 - 8.0f / 15.0f);
        return sum;
    }

    /**
    * h function in Pdm nonlinear material is \f$ \sum^4_{i=1}{p_iA_i} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T hFunction(T x, T* ptr_material)
    {
        return gFunction(x, ptr_material);
    }

    /**
    * first derivative of g function in Pdm nonlinear material is 2*(x-1)
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T dgFunction(T x, T* ptr_material)
    {
        T nonlinear_params[4] = { ptr_material[2],ptr_material[3] ,ptr_material[4] ,ptr_material[5] };

        auto x2 = x * x;
        auto x3 = x2 * x;
        auto x4 = x2 * x2;

        auto sum = nonlinear_params[0] * (x - 1.0f / x);
        sum += 0.5f*nonlinear_params[1] * (x2 - 1.0f / x2);
        sum += 1.0f / 3.0f*nonlinear_params[2] * (x3 - 1.0f / x3);
        sum += 0.25f*nonlinear_params[3] * (x4 - 1.0f / x4);
        return sum;
    }

    /**
    * first derivative of h function in Pdm nonlinear material is 2*(x-1)
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T dhFunction(T x, T* ptr_material)
    {
        return dgFunction(x, ptr_material);
    }

    /**
    * second derivative of g function in Pdm nonlinear material
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T ddgFunction(T x, T* ptr_material)
    {
        T nonlinear_params[4] = { ptr_material[2],ptr_material[3] ,ptr_material[4] ,ptr_material[5] };

        auto x2 = x * x;
        auto x3 = x2 * x;
        auto x4 = x2 * x2;
        auto x5 = x4 * x;
        auto sum = nonlinear_params[0] * (1.0f + 1.0f / x2);
        sum += nonlinear_params[1] * (x + 1.0f / x3);
        sum += nonlinear_params[2] * (x2 + 1.0f / x4);
        sum += nonlinear_params[3] * (x3 + 1.0f / x5);

        return sum;
    }

    /**
    * second derivative of h function in Pdm nonlinear material
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T ddhFunction(T x, T* ptr_material)
    {
        return ddgFunction(x, ptr_material);
    }

    /**
    * anisotropic function for Pdm nonlinear material, it is isotropic.
    * @param[in] dir  bond direction in material space
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T anisoWeight(glm::tvec3<T> dir, T* ptr_material)
    {
        return 1.0f;
    }
};

}